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We evaluate analytically and numerically the size of the frozen core and various scaling laws for 
critical Boolean networks that have a power- law in- and/or out-degree distribution. To this purpose, 
we generalize an efficient method that has previously been used for conventional random Boolean 
networks and for networks with power-law in-degree distributions. With this generalization, we can 
also deal with power-law out-degree distributions. When the power-law exponent is between 2 and 3, 
the second moment of the distribution diverges with network size, and the scaling exponent of the 
nonfrozen nodes depends on the degree distribution exponent. Furthermore, the exponent depends 
also on the dependence of the cutoff of the degree distribution on the system size. Altogether, we 
obtain an impressive number of different scaling laws depending on the type of cutoff as well as 
on the exponents of the in- and out-degree distributions. We confirm our scaling arguments and 
analytical considerations by numerical investigations. 



I. INTRODUCTION 

Boolean networks are often used as generic models for 
the dynamics of complex systems such as social and eco- 
nomic networks, neural networks, and gene or protein 
interaction networks [1, 2]. Whenever the states of the 
constituents of the system can be reduced to being ei- 
ther "on" or "off" without loss of important information, 
a Boolean approximation captures many features of the 
dynamics of real networks [3]. In order to understand 
the generic behavior of Boolean networks, random mod- 
els were investigated in depth [4], although it is clear 
that neither the connection pattern, nor the usage of 
update functions of biological networks is reflected re- 
alistically in such random models. For random models, 
the distinction between frozen, chaotic, and critical net- 
works has become commonplace [5, 6], with frozen net- 
works having short attractors where almost all nodes are 
fixed at one value, while chaotic networks have attrac- 
tors the length of which increases exponentially with the 
system size. Critical networks, which are "at the edge 
of chaos", are considered particularly relevant for bio- 
logical systems, and for this reason many investigations 
have concentrated on such critical networks. One impor- 
tant result of these investigations was that the number 
of nodes that do not become frozen on the attractors 
increases as TV 2 / 3 with the network size TV [7-9]. 

The majority of recent research on Boolean networks 
was devoted to networks with more realistic features. In 
particular, many natural networks are scale free, which 
means that the number of connections per node follows 
a power law [10]. Typically, the power-law exponent of 
the degree distribution is between 2 and 3, which means 
that the second moment of this distribution diverges with 
network size. Due to their relevance to natural systems, 
scale-free Boolean networks have been investigated by 
various authors. In order to keep these models as sim- 
ple as possible, connections between nodes are made at 
random within the constraints given by the degree distri- 
bution. The Boolean functions are usually also assigned 



at random from a certain set of functions, e,g., thresh- 
old functions. Observations made in computer simula- 
tions for these networks are that attractors are shorter 
and frozen nodes are more numerous in critical scale- 
free networks compared to conventional random Boolean 
networks with the same total number of links and of 
nodes [11, 12], that attractors are sensitive to pertur- 
bations of highly connected nodes, but not of sparsely 
connected nodes [12, 13], and that scale-free Boolean 
networks evolve much faster and more steadily towards 
a target pattern of an "output" node than conventional 
random Boolean networks [14]. Analytical results are 
mostly limited to calculating the phase diagram using 
the annealed approximation [13, 15, 16]. Similarly to con- 
ventional random Boolean networks, which have a fixed 
in-degree, scale- free Boolean networks show also frozen, 
critical, or chaotic dynamics, depending on the parameter 
values. A comparison between results obtained from the 
annealed approximation and from computer simulations 
is performed in [17], giving not always an agreement be- 
tween the two approaches. Three years ago, Drossel and 
Greil [18] derived analytically the scaling exponents for 
the number of nonfrozen nodes in critical Boolean net- 
works with a scale-free in-degree distribution. The values 
of these exponents depend continuously on the exponent 
of the in-degree distribution, if it is smaller than 3. So far, 
the equivalent case of scale-free out-degree distributions 
was not yet investigated, although it is relevant for natu- 
ral systems [10]. From computer simulations of Boolean 
networks with a scale- free out-degree distribution, it is 
known that the properties of attractors are different from 
the case of a scale-free in-degree distribution [19]. 

It is the purpose of this paper to derive the scaling 
exponents of the number of nonfrozen nodes for Boolean 
networks that have a power-law out-degree distribution. 
To this goal, we generalize an efficient method that has 
previously been used for conventional random Boolean 
networks [8, 9] and for networks with power-law in-degree 
distributions [18]. The result is a stunning variety of scal- 
ing laws depending on the in- and out-degree exponent as 
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well as on the type of cutoff used in both cases. In partic- 
ular, we also find that the dependence of the scaling ex- 
ponents on the degree distribution shows opposite trends 
for scale-free in-degree and for scale-free out-degree dis- 
tributions. For one of the many cases investigated by 
us, we obtain the scaling law given in [20]. However, we 
find this scaling law for a different situation than the one 
considered in [20]. We confirm our analytical results by 
computer simulations. 



II. MODEL 

We consider Boolean networks consisting of TV nodes, 
where each node i has a Boolean value g\ G {0, 1} and is 
connected to other nodes via its in- and outputs. Further- 
more, each node is assigned a Boolean update function. 
Connections and update functions are chosen at random, 
given the distribution P(k- m ) and P(k out ) of the number 
of inputs and outputs of the nodes, and the distribution 
of update function. 

In conventional random Boolean networks, all nodes 
have the same number k of inputs, P(k{ n ) = Sk iT1 ,k and a 
Poisson distribution of the number of outputs, 

Z^&out 

P(k oat ) = —e- k , (1) 

since each node receives its input from each other node 
with the same probability. Different probability distribu- 
tion for Boolean functions are used, with biased functions 
being a familiar choice, where for each possible combina- 
tion of the ki input values the output is 1 with probability 
p and with probability 1 — p. By adjusting the value 
of p, the network can be made critical. In this paper, 
we use only constant and reversible functions. For each 
value of the in-degree, there are two constant functions, 
which take the same value for all possible inputs, and two 
reversible functions, which are defined by the condition 
that the change of one input always changes the output. 
For k [n = 2, these functions are XOR and NOT XOR. 
We will argue in the end that our results are also valid 
for other choices of update functions, in particular for 
biased functions. 

In this paper, we consider scale-free networks, where 
either the distribution of inputs is a power law, P(h n ) ex 
^in 7in or the distribution of outputs is a power law, 
P(k out ) oc k~^ t ou \ or both are a power law. We only 
consider the case that there is no correlation between 
the in-degree and the out-degree of a node. We focus on 
the interesting case of 7i n , 7 ou t G (2, 3), where the second 
moment of k[ n or fc out diverges, but the mean value is 
well defined. The mean values of the in- and out-degree 
distribution have to be identical, since the total number 
of inputs and outputs must be the same because each 
link connects an input with an output. 

We generated scale-free in- or out-degree distributions 
P(k) oc /c -7 in two ways, which lead to a different scaling 



of the cutoff values & max (7V) with network size N. One 
the one hand, we draw each value ki at random from the 
distribution P (fc), which leads to a cutoff 

fc ffiax ex N . (2) 

Alternatively, we fixed the number of nodes with the in- 
or output value k to exactly c • N • P (k) rounded to the 
next integer, while adjusting c ^ 1 such that the size of 
the sample is as close as possible to N. This leads to a 
cutoff 

P x ex N^f . (3) 

We fixed the minimum value of k in scale-free distribu- 
tions to 2, since this leads to more nodes with larger 
values of k and therefore to a faster approach to the 
asymptotic behavior with increasing N. Furthermore, 
we adjusted the number of nodes with k = 2 connec- 
tions such that the mean value of k does not change with 
N, but equals the asymptotic value taken in the limit 
N — >> oo. This also reduces finite-size effects. 

In order to make the networks critical, the proportion 
r c of reversible functions was chosen such that the change 
of the state of one node propagates on average to one 
other node, implying (1 — r c ) -0 + r c • (k) — 1, and leading 
to 

For such critical networks, the number of nodes that do 
not become frozen when the system is on an attractor 
increases sublinearly with network size, with a power-law 
exponent that will be determined below. 



III. ALGORITHM: DETERMINING THE 
FROZEN CORE STARTING FROM CONSTANT 
FUNCTIONS 

An elegant way to determine the frozen core was sug- 
gested in [8, 9] and is generalized in the present paper 
such that it can be used for out-degree distributions that 
are not Poissonian. This method is based on the as- 
sumption that almost all frozen nodes can be obtained 
by starting from the nodes with a constant update func- 
tion and by determining iteratively all nodes that become 
frozen because some of their inputs are frozen. This as- 
sumption is valid in many cases, in particular for net- 
works that have only constant and reversible functions. 
If the assumption is not correct, one can nevertheless 
expect the same scaling exponents, but must use other 
methods to obtain the frozen core [21]. The main idea 
of the method is to not specify the network in advance, 
but to choose the connections within the network while 
determining the frozen core. We used the algorithm for 
our analytical calculations as well as for our computer 
simulations. The steps of the algorithm are as follows: 
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1. Initially, each of the TV nodes is assigned an up- 
date function, a number of inputs, and a num- 
ber of outputs according to the rules given above. 
These three assignments are made independently 
from each other, in order to avoid correlations be- 
tween them. Nodes with a reversible function are 
placed into "containers" d according to their num- 
ber ki n = i of inputs. Nodes with constant func- 
tions are put into the container Cq. We denote the 
number of nodes in container d by \d\. Further- 
more, we denote the total number of nodes in all 
containers by Nf = X^=cT l^l> anc ^ the total num- 
ber of inputs to nodes that are not in containers d 
with i > 1 by kf n . The initial value of kf n is identi- 
cal to the number of inputs to nodes with constant 
functions. All these values will change during the 
algorithm. In particular, iVf will decrease by 1 with 
each step. 

2. Then, the following steps are iterated until |Co| = 
0: 

(a) Select one node from container Co- Its number 
of outputs is denoted by k^ t . 

(b) Choose at random k^ t out of all ^ i \Ci\-\-kf n 
inputs to become connected to the outputs of 
the selected node. 

(c) If m > inputs of a node in a container i > 1 
became connected to the selected node, move 
this node from its container d to d-m- 

(d) Reduce the number kf n of unconnected inputs 
to nodes with constant functions by the num- 
ber of those that became connected to the se- 
lected node in Co- This ensures that the to- 
tal number of outputs in all containers always 
equals the total number of remaining inputs, 

(e) Remove the selected node from the system. 
This implies the replacement iVf := iVf — 1. 

3. The final value of iVf is identical to the number 
of nodes that do not belong to the frozen core of 
the particular network that was constructed by per- 
forming this algorithm. The probability distribu- 
tion of iVf follows from the stochastic process im- 
plemented in this algorithm. 

In order to finish the construction of the network, the 
inputs and outputs of the remaining nodes, which con- 
stitute the nonfrozen part of the network, should also be 
connected at random. However, since we are only inter- 
ested in the number of nonfrozen nodes and not in the 
attractors of the networks or the structure of relevant 
components, we omit this step. 



IV. ANALYTICAL CONSIDERATIONS 

Based on the algorithm outlined in the previous sec- 
tion, the scaling of the number of nonfrozen nodes with 
network size can be determined analytically. To this pur- 
pose, we define the parameter e = which is the pro- 
portion of nodes that have not yet become frozen. Ne- 
glecting fluctuations, the number of nodes in containers 
i > 1 can be expressed as [18] 

* max / 7 \ 

\C i \ = J2\Cr\e i (l-e) l -^.j , (5) 

for i > 1, where \Cj m | is the number of nodes in container 
I at the end of step 1. of the algorithm. When e is small, 
only a small proportion of all inputs are still present, 
and most nodes that are in container i were initially in 
containers with values I ^> i. Therefore, expression (5) 
can be approximated for small e by 

/*max 
\C'f\e- d rdl. (6) 

Due to the condition Eq. (4) for critical networks, the 
initial number of nodes in container Co is 

\C™\ = ^£(i-l)\Cf \ , (7) 

which is equivalent to the condition 

N t = ^2i\d\, (8) 

i=0 

which holds during the entire algorithm (again neglecting 
fluctuations), since on average one input in the containers 
i > 1 will become connected to the selected node during 
one step, because each of the Nf(k) remaining inputs 
connects to the (k) outputs of the selected node with the 
same probability. Consequently we have 

|Co|=5>-l)|Q| (9) 

i=l 

during the entire algorithm, and the contents of all con- 
tainers will reach the value zero at the moment when e 
reaches zero. 

For small values of e, the sum in Eq. (9) is dominated 
by the i = 2 term, as can be concluded from Eq. (5). 
This means that for small e almost all non-frozen nodes 
are in C\ and C2 and that 

\Co\ * \C 2 \ (10) 

apart from fluctuations. 

Due to random fluctuations, |Co| will typically reach 
the value zero while Nf is still larger than zero. This will 
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happen when the standard deviation of |C 1 , denoted as 
<7|C |j becomes comparable to the average value of |Co|, 
which in turn is identical to the average value of | C2 1 5 
giving the condition for the end of the process 

IOjI^ici- (ii) 

In order to obtain from this relation a condition for the 
scaling of the final value Nf with TV, we must express 
both sides in terms of N and Nf (or, equivalently, e). 

In the following, we determine the variance 0\c Q \ °f 
nodes in container Cq. In every step exactly one node 
is removed from the system and some of the nodes from 
Ci with i > are moved to Cq. The total number of 
nodes in all containers, Nf + |Cb| decreases exactly by 1 
for each iteration of the algorithm. Fluctuations of the 
total number Nf of nodes in containers Ci with i > 
are deviations from the mean values given by (6). These 
deviations are identical, but with opposite sign, to fluctu- 
ations of the total number of nodes in container Co, since 
the total number of nodes in the containers decreases in 
a deterministic manner and can therefore not show ran- 
dom fluctuations. Since for small e the vast majority of 
remaining nodes are in Ci, the fluctuation of the num- 
ber of nodes in containers Ci with i > is dominated by 
those in container C±, leading to 

°"|C | = (J N i ^ CT|d|. (12) 

There are two contributions to this variance: 

1. The number of nodes Nf remaining in containers 
Ci with i > 1 has to be on average the number of 
remaining outputs divided by their mean degree. 
This gives the following contribution to the vari- 
ance 

2. The second contribution comes from the fact that 
for small e the number of nonfrozen inputs is Pois- 
son distributed with a variance identical to the 
mean value, which is proportional to Nf. Since the 
vast majority of nodes only have one input, this is 
also the variance of \C\\, leading to a 2 ^ oc Nf. 

For networks with a scale-free out-degree distribution 
with an exponent between 2 and 3, the first term dom- 
inates and is proportional to Nf multiplied by a power 
of N . Otherwise, the first and second term give together 
cr^ f oc Nf. We conclude that 

(T\c Q \ oc VNf<T kottt oc ^NfN A (14) 

with an exponent A that depends on the out-degree dis- 
tribution. If the variance of the out-degree distribution 
is finite, we have A = 0. This holds when the out-degree 
of all nodes is identical or is Poisson distributed, or when 
it is a power-law distribution with an exponent 7 ou t > 3. 



When the out-degree distribution is a power law with 
7out G (2,3), we have 

<x (Cut x ) 3_7out , (15) 

leading with Eq. (2) or Eq. (3) to A = 3 - 7 ou t or A = 
(3 — Tout) /Tout- For Tout = 3, all three expressions agree 
with each other and give A = 0, as it must be. 

In order to complete the calculation, we need an ex- 
pression for I C 2 1 in Eq. (11). This expression is contained 
in Eq. (6) and was already given in [18] for the case of a 
Poissonian out-degree distribution. When the in-degree 
distribution has a finite second moment, which is the case 
for fixed, Poissonian, or power-law in-degree distributions 
with Tin > 3, the integral in Eq. (6) converges even when 
the cutoff is set to infinity and e is set to 0, leading to 

|C 2 | oc TVe 2 . (16) 

For networks with a scale-free in-degree distribution we 
can rewrite Eq. (6) as 

\d\ ex e l N ( ^ e- el l l -^dl . (17) 

J i 

In the case Tin > 3, the integral converges, leading again 
to Eq. (16). In the case Tin £ (2,3), the second moment 
diverges, and the value of the integral is determined ei- 
ther by the upper limit of the integral fc£J ax or by the in- 
verse exponential decay constant 1, whichever is smaller. 
This gives 

\C 2 \ oc Ne 2 (min (k™*, ^ ^ . (18) 

For the case ^ ax > 1, we obtain 

|C 2 | ociVe^- 1 , (19) 

independently of the scaling of /c[JJ ax with N. For /cJJJ ax < 
1 , we obtain for /c[JJ ax oc TV 

\C 2 \ oc7V 4 -^e 2 , (20) 
and for A:^ ax ex TV 1 / 7 - 

\C 2 \ ocTvie 2 . (21) 

Fo r Tin = 3, all four expressions for \C 2 \ give the same 
result, as it must be. 

Taking the four cases together, we can write 

\C 2 \ (xN a e b (22) 

with different values for a and b. 

Inserting this result together with Eq. (14) into 
Eq. (11), we obtain the desired scaling law of the num- 
ber of non-frozen nodes 7V n f, which is identical to the 
final value of Nf, as 

7V nf cx eN cx N 1 = N x . (23) 
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Since both /c[JJ ax and | scale with N in a nontrivial way, 
and since the scaling of e with N depends on the result 
of Eq. (18), it is not always possible to tell in advance 
whether /c[JJ ax is larger or smaller than K In such cases, 
we assumed first one version of the inequality, and if this 
gave an inconsistency, we used the other version. If we 
write /c[JJ ax ex N z and use the fact that e ex N x ~ l (with 
x defined in eq. (23)), the consistency condition reads 



C ax ^ - z 



1^0. 



(24) 



In the standard case, where the second moments of the 
in- and out-degree distributions are finite, we have A = 0, 
a = 1, and 6 = 2, leading to the well-known result 



nf OC 



AT 3 . 



(25) 



Since we consider three cases for both the in-degree and 
the out-degree distribution, there are altogether 9 differ- 
ent relations for the scaling of the number of non-frozen 
nodes with the total number of nodes in a critical net- 
work. These 9 cases are listed in Tab. I and represented 
graphically in Fig. 1. The most striking feature of these 
results is that the scaling exponent for the number of non- 
frozen nodes increases with 7i n , but decreases with 7 ou t- 
In the case of a scale-free out-degree disribution, the ex- 
ponent characterizing the proportion of nonfrozen nodes 
that have two nonfrozen inputs increases with decreasing 

Tout- 

There are several ways to check the correctness of these 
expressions: 

• Every expression for x must give | for 7i n = 7 out = 
3. If both distributions are scale free, x must take 
the expression for the corresponding Poisson distri- 
bution when one of the 7 values is set to 3. 

• x must be in the interval [0, 1] for 7 G [2, 3]. 

• For a scale-free in-degree distribution, one of the 
two possible expressions for x should fulfill the con- 
sistency condition Eq. (24). If both fulfill the con- 
dition, the values of x must be identical in the two 
cases. 

• From Eq. (24) follows that x = 1 — z fo the border 
case. For z = 1 (/c£J ax cx N), this meas that x = 
is only possible on the border of the considered 



range. For z 



the values of x must be x 



1 — — on the border, independently of the out- 
degree distribution. 

We also did the ultimate check, which is computer simu- 
lations, as described in section VI. 

In a similar way, one obtains scaling laws for | C2 1 - To 
this purpose, the scaling of e with N from Eq. (23) must 
be combined with Eq. (22), giving 



|C 2 | ocTV 



a+b- 
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FIG. 1. (Color online) Graphical representation of the nine 
cases listed in Table 1. Top: exponent x. Bottom: exponent 



The values of y for the nine cases are also listed in Tab. I, 
and I is visualized in Fig. 1. The check for correctness 
can be made in the same way as for the exponent x, with 



(26) the only difference that y must be | for 7 = 3. 
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TABLE I. The 9 different expressions for the exponents x and y that characterize the scaling N n { oc N x of the number of 
nonfrozen nodes N n f with N and the scaling of the number of nonfrozen nodes with two nonfrozen inputs |C 2 | oc N y with 
N. The expressions obtained in the special case 7m = 7out are also specified. Where necessary, two expressions for x and the 
corresponding conditions are given, as well as the boundary case for which both expressions hold simultaneously. 



V. EXTENSION TO OTHER SETS OF UPDATE 
FUNCTIONS 

All results so far were derived by using only constant 
and reversible update functions. However, the algorithm 
can be generalized to more general sets of update func- 
tions, in particular to biased functions. For general sets 
of functions, in step 1. nodes with constant functions are 
placed in container Co, while all other nodes are placed 
in containers according to their number of inputs. Also 
step 2. has to be modified. When m > 1 inputs of a node 
in container i become connected to the selected node, 
the node in container i may freeze completely and is then 
moved to container Co instead of Ci- m . This occurs with 
a probability that depends on the set of update functions 
and is given by the probability that a randomly chosen 
function of i inputs becomes constant when m of the in- 
puts become frozen. In this case also the value of fcP n has 
to be increased by i — m. Clearly, the container method 
only works when the probability distribution of the up- 
date functions with i inputs is identical to the conditional 
probability distribution that is obtained by freezing I — i 
inputs of nodes that are initially in a container C\ with 
I > i. 

The analytical considerations become slightly different, 
but lead to the same conclusions. Eq. (5) still states that 
towards the end of the algorithm only containers Co, Ci, 
and C2 need to be considered. Nodes in Co or C\ will 
behave in the same way as for the case of constant and 



reversible update functions. The only difference is that 
a certain proportion of the nodes in C2 become already 
frozen if only one of their inputs is connected to a frozen 
node. Such nodes could in principle be placed in con- 
tainer Ci, and then the situation would be identical to 
the one of only constant and reversible functions. Due 
to the fact that \C\\ ^> | O2 1 , this make no significant dif- 
ference to the calculations and therefore has no effect to 
the scaling laws. 

VI. COMPUTER SIMULATIONS 

In order to confirm our analytical calculations, we per- 
formed computer simulations of the algorithm. Instead 
of a Poisson distribution we used a fixed value for k. If 
the value of (k) required for criticality was not an integer, 
we used a mixture of the two neighbouring integer val- 
ues for generating the distribution of k values. Avoiding 
unwanted correlations between in-degree and out-degree 
distribution or between a degree distribution and the dis- 
tribution of Boolean functions in the case of scale-free dis- 
tributions turned out to be challenging. The simulations 
were done in the following way. For each set of in- and 
out-degree distributions and each value of 7, we run the 
algorithm between 10 3 and 10 6 times (the smaller value 
applies to larger values of N, due computation costs) for 
values of N e {2 10 , 2 11 , . . . , 2 20 , . . .} and thus obtained a 
distribution for the number of nonfrozen iV n f and the final 
size of nodes in container C2. The data were smoothened 
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FIG. 2. Typical example for the probability distribution of 
nonfrozen nodes. This sample is for critical networks with a 

scale free in- and out-degree distribution with 7 = 2.3 and 

1 

a cutoff scaling as /c max oc N 7 . The number of realisations 
was 2 • 10 5 , and number of nodes in the network N = 2 10 for 
the solid and N = 2 11 for the dotted curve. The curves were 
collapsed using the exponent 0.635, which is the best value 
our algorithm found for these two curves. 




FIG. 3. (Color online) Example for the color coding in Fig. 4. 
Each subplot is similar to Fig. 2. The only difference is the 
number of nodes, which was in this case 2 17 and 2 18 . The 
numbers above the graphs give the exponents used for scal- 
ing the two curves, and the color indicates the quality of the 
collapse. The theoretical exponent is used in the last plot. 



by performing a logarithmical binning with a step width 
of 1.05. To further smoothen the first part of the curve 
where the slope is very small (in the log-log plot), we 
combined groups of neighboring bins that comprise ap- 
proximately the same number of events. A typical result 
is shown in Fig. 2, where the axes were scaled using the 
exponent that gives the best data collapse. 

However, the quality of the data collapse is not always 
that good, in particular when 7 is close to 2 or when 
N is small. In order to compare the analytical expres- 
sion, which should become exact in the limit N — » 00, 
to the simulation results, we established an automated 
procedure for quantifying the quality of the date col- 
lapse between two curves for system sizes Ni = 2 n and 
N2 = 2 n+1 for the entire possible interval (0, 1). This au- 
tomated procedure gives an optimal (TV-dependent) value 
of the exponent, which should approach the theoretical 
value as N increases. The quality of the collapse was 
quantified by a fitness value, which is the mean of the 
absolute value of the distance between the two linearly 
interpolated curves, omitting the first 20 percent of the 
curves. We convinced ourselves that higher fitness val- 
ues correspond to what is intuitively considered a better 
collapse. The results are shown in Fig. 4, using colors to 
indicate the fitness. Each box shows on the y axis the 
scaling exponent used for the collapse and on the x axis 
the value of Ni of the pair that was compared. The color 
scale was set such that the interval between the mini- 
mum and maximum fitness value was stretched to [0, 1] 
and divided by 0.7 and then taken to the power of 0.5 



to improve color resolution for minimal (optimal) fitness 
values close to 0. Fig. 3 gives an impression of the quality 
of the collapse associated with the different colors. Some 
of the boxes contain white areas due to missing simula- 
tion values or to large fluctuations due to poor statistics. 
The theoretical value for the exponent lies in almost all 
boxes in the area of best fitness values. In some of the 
cases one can see that this area is moving down or up 
with increasing system size N±. This indicates finite size 
effects. 

We also did this evaluation for the scaling exponents 
for I C 2 1 (not shown), and we found again that the agree- 
ment between theory and computer simulations is very 
good. 

VII. CONCLUSIONS 

Using analytical calculations and computer simula- 
tions, we have determined the scaling of the number 
of nonfrozen nodes of critical scale-free random Boolean 
networks with network size, leading to a stunning variety 
of different scaling laws. Our calculations are based on an 
algorithm that determines the frozen core of the network 
by an iterative procedure. The scaling exponent for the 
number of nonfrozen nodes does not only depend on the 
values of two exponents of the power-law degree distri- 
butions (if they are smaller than 3) , but also on whether 
these scale-free distributions are generated by randomly 
sampling the degree of each node from a power-law dis- 
tribution or by exactly matching the degree distribution 
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of each realization of the network to the desired power 
law. We calculated the scaling laws by generalizing a phe- 
nomenological theory that has been used earlier for criti- 
cal random Boolean networks with Poissonian out-degree 
distributions. Furthermore, we performed computer sim- 
ulations using a very efficient algorithm that allowed us 
to test the theoretical results for network sizes larger than 
2 20 . These computer simulations confirm our analytical 
considerations. This work thus fills an important gap in 
our understanding of scale-free critical random Boolean 
network. 

Our results show that the size of the nonfrozen part of 
the critical network decreases with decreasing 7i n but in- 
creases with decreasing 7 ou t- I n the case where the cutoff 
of the out-degree distribution scales as TV, the scaling ex- 
ponent of the nonfrozen part of the network approaches 1 
as 7 out decreases towards 2. This means that a finite pro- 
portion of all nodes remain nonfrozen in this limit case. 
In the case 7i n = 7 out = 7, the trend of the scaling expo- 
nents with 7 depends on the scaling of the cutoffs with 
network size. The opposite trends observed for scale- free 
in- and out-degree distributions can be explained as fol- 
lows: a smaller value of 7i n leads to smaller fluctuations 
in the contents of the containers (i.e., in the growth of the 
frozen core) and therefore to a later stop of the freezing 
algorithm, while a smaller value of 7 ou t leads to larger 
fluctuations in the contents of the containers and there- 
fore to an earlier stop of the freezing algorithm. 

The proportion of nonfrozen nodes with two nonfrozen 
inputs increases as 7 ou t decreases towards 2, and it ap- 



proaches the value 1 in the case where the cutoff of the 
out-degree distribution is proportional to N. In contrast, 
in the case of Poissonian out-degree distribution the num- 
ber of nonfrozen nodes with two nonfrozen inputs scales 
as the square root of the total number of nonfrozen nodes, 
irrespective of the in-degree distribution. For this case, it 
was shown in [22] that the computational core of the net- 
work (i.e., the set of relevant nodes), which determines 
the number and length of attractors, scales also as the 
square root of the number of nonfrozen nodes. Our re- 
sults show that the computational core increases when 
the out-degree distribution becomes scale free, and the 
majority of relevant components are not any more sim- 
ple loops. This means that the length of attractors is 
much larger for scale-free out-degree distributions than 
for scale-free in-degree distributions. 

Finally, we want to emphasize that we only investi- 
gated situations where the in- and out-degree of a nodes 
were uncorrelated. However, the case where the two de- 
grees are correlated (they can for instance be identical) 
is also relevant. It applies in particular to networks that 
have undirected connections, which means that the in- 
and outgoing connections are identical. As argued in 
[16], the undirected case with a degree-distribution ex- 
ponent 7 corresponds to the uncorrelated directed case 
with an exponent 7 — 1. 
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FIG. 4. (Color online) The quality of the data collapse (color coded) of the distributions P(N n f) obtained for two N values 
that differ by a factor of 2, as a function of N (x-axis log. 10 3 to 10 8 ) and the exponent used for the collapse (y axis lin. 
to 1). The darker the color, the better the collapse. The different plots are for different combinations of the in- and out- 
degree distributions. When the in-degree distribution was identical to the out-degree distribution, we generated the out-degree 
distribution by using exactly the same sample as the one obtained for the in-degree distribution. The lines give the value 2/3 
(dotted line), the theoretically predicted exponent (dash-dotted line), and the best exponent for the collapse of the pairs of 
curves (solid line). 



